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ABSTRACT 

We present an exact three-dimensional wave solution to the shearing sheet 
equations of motion. The existence of this solution argues against transient 
amplification as a route to turbulence in unmagnetized disks. Moreover, because 
the solution covers an extensive dynamical range in wavenumber space, it is an 
excellent test of the dissipative properties of numerical codes. 

Subject headings: accretion, accretion disks; instabilities; solar system: formation 

1. Introduction 

The combination of magnetic fields and differential rotation is a volatile one, leading to 
the magnetorotational instability, or MRI, under very general conditions (e.g. Balbus 2003). 
Laboratory evidence for the MRI in liquid sodium has recently been claimed (Sisan et al. 
2004), and there can be little doubt that in magnetized disks the MRI will be present and a 
source of the enhanced angular momentum transport needed for accretion to proceed. But 
not all astrophysical disks need be everywhere ionized at a level at which the field and the 
gas are well-coupled. Therefore the question of whether or not an unmagnetized, laminar gas 
in Keplerian rotation is vulnerable to turbulent breakdown continues to attract attention. 

The nature of the local stability of Keplerian disks has been investigated by detailed 
numerical simulation (Balbus, Hawley, & Stone 1996; Hawley, Balbus, & Winters 1999). 
There is no evidence of nonlinear instability. Indeed, the notion that Keplerian disks are 
intrinsically turbulent would be most surprising and problematic to MHD numericists, who 
routinely check the validity of an MHD code by ensuring that the disk quickly returns to 
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stability when the magnetic field is removed. Even when nonlinear perturbations are directly 
seeded, they do not lead to turbulence. 

Numerical simulations have been criticized, however, on the grounds that they lack 
the necessary resolution to reveal nonlinear instabilities (e.g. Richard & Zahn 1999; Af- 
shordi, Mukhopadhyay, & Narayan 2005). The claim is that the effective Reynolds number 
of available codes is too low, and that the numerical schemes dissipate nascent nonlinear 
instabilities. In some cases, a specific destabilizing mechanism is suggested. Afshordi et al. 
(2005) suggest, for example, that the transient growth exhibited by linear nonaxisymmetric 
perturbations, in the course of evolving from leading to trailing wavenumbers, might become 
sufficiently large to cross a threshold and initiate nonlinear turbulence. 

The first systematic numerical studies of Keplerian stability are now a decade old, and 
numerical codes have developed and improved. Our physical understanding of disks has also 
deepened. We are motivated to reexamine the problem of the stability of Keplerian disks. 
In this paper, we present precise analytic solutions for three-dimensional incompressible per- 
turbations in Keplerian disks, recovering the previously known two-dimensional solutions as 
special cases. We show that in fact all such solutions are valid both linearly and nonlinearly 
because the nonlinear terms in the equations of motion vanish. Thus, these solutions by 
themselves do not trigger nonlinear instabilities, regardless of the level of transient amplifi- 
cation attained. Since the solutions exhibit both strongly leading and strongly trailing large 
wavenumbers, they cover a very large dynamical range. They are, in fact, a good test of 
the dissipation properties of numerical codes. Conversely, the numerical codes test the sta- 
bility of these exact time-dependent solutions, which would otherwise be a difficult analytic 
calculation. 

Recently, laboratory experiments have weighed in on the controversy. A high Reynolds 
number (~ 10 6 ) Couette flow experiment found that Keplerian-like profiles are nonlinearly 
stable (Ji 2006, private communication). This finding is compatible with numerical simu- 
lations, but in conflict with an earlier experiment (Richard 2001) that claimed that such 
profiles were unstable. The Ji experiment uses differentially rotating endcaps to break up 
Ekman boundaries (which otherwise cause unwanted interior circulation), and Doppler ve- 
locimetry techniques to verify that the rotation profiles matched their analytic counterparts. 
The Reynolds stress in the Ji experiment has also been measured, and residual fluctuations 
in the radial and azimuthal velocities are uncorrelated: the tensor vanishes on average. 

These results are compelling reasons to believe that local Keplerian flow is nonlinearly 
stable, and that the numerical simulations have not been prohibitively dissipative or oth- 
erwise qualitatively misrepresentative. The explicit solutions we present in this paper will 
provide numericists with a well-posed numerical test bed that will sharpen our understanding 
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of local small scale behavior in Keplerian disks. 

An outline of this paper is as follows. In §2 we provide an overview of the issues that 
have been debated. Section 3 sets up the mathematical background, and the full analysis is 
given in §4. The solutions derived in §4 are applied as a code test problem in §5. Finally, §6 
summarizes our work. 



Axisymmetric waves propagating in the interior of an unmagnetized disk are stable if 
and only if the specific angular momentum / increases with cylindrical radius R. This is the 
classical Rayleigh stability criterion, and it is generally satisfied in astrophysical disks. When 
stable, fluid displacements in the disk plane oscillate about their equilibrium circular orbits 
with a characteristic frequency k. The value of k 2 is directly proportional to the specific 
angular momentum gradient I. In particular, 



where R is the cylindrical radius and il = l/R 2 is the local disk angular velocity. The formal 
transition to a negative value of n 2 corresponds to the onset of rotational instability. It is to be 
noted that this result applies only to axisymmetric disturbances; there is at present no known 
rigorous nonaxisymmetric criterion. (Hence the ongoing debate.) A radially discontinuous 
rotation profile, for example, would be unstable to non-axisymmetric disturbances by the 
Kelvin-Helmholtz instability, even if the rotation rate increased outward and was formally 
Rayleigh-stable. 

For a Keplerian disk, k 2 = Q 2 , which is a stable profile by the Rayleigh criterion. 
Absent pressure forces, this is unconditionally stable: the behavior of Keplerian orbits is 
not an object of current debate. Pressure forces must, therefore, be at the heart of any 
proposed instability mechanism. This is already a potential difficulty for proponents of 
instability, since any system that is close to Keplerian must be characterized by very small 
radial background pressure gradients, and perturbative pressure gradients generally behave 
acoustically. Even global instabilities in flows that would be locally Rayleigh-stable revert to 
stable behavior as the disk approaches a Keplerian profile (see, e.g., Goldreich, Goodman, 
& Narayan 1986). 

The claims to nonlinear instability discussed in the Introduction are of course meant to 
apply to systems which are linearly stable. It is not unprecedented for linear waves to be sub- 
ject to a finite amplitude instability (e.g., Infeld & Rowlands 1990), and even at arbitrarily 
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small amplitudes, one-dimensional adiabatic sound waves show nonlinear steepening behav- 
ior. (In fact, the shearing wave solutions to be discussed are, in their tightly-wrapped phase, 
subject to dissipative instabilities that may well prevent them from achieving even modest 
amplitudes.) But proponents of nonlinear Keplerian instability are arguing for something far 
more novel than wave modification. The claim is that well-behaved local linear disturbances 
that formally begin and end their existence as small amplitude waves, heavily dominated by 
restoring forces, in actuality become self-sustaining turbulence. If only the small nonlinear 
terms were properly resolved, it is argued, the simulations would reveal this behavior. In 
this view, current disk codes would be not merely imprecise, but grossly misrepresentative 
of the most fundamental character of Keplerian flow. 

Planar Couette and Poiseuille viscous flows are two classic shear profiles that exhibit a 
breakdown into turbulence. Their behavior has motivated the notion that Keplerian rotation 
ought to do the same. Planar Couette flow is linearly stable (in a marginal sense) but 
nonlinearly unstable, whereas Poiseuille flow is linearly unstable if the appropriately defined 
Reynolds number Re exceeds a large critical value (Re cr u = 5772, see e.g., Drazin & Reid 
[1981] for details). Both flows are exquisitely sensitive to nonlinear disturbances. 

Planar Poiseuille and Couette flows have boundary layers and linear dynamical responses 
lacking a dominant restoring force. The presence of viscous boundary layers means that the 
Reynolds number is an important feature of the flow stability, and at the same time allows for 
linear instability in Poiseuille flow at sufficiently high Reynolds number. By way of contrast, 
local Keplerian flow has no viscous boundary layers, and its destabilizing shear is strongly 
stabilized by the Coriolis force, 2Q > \dQ/d\nR\. It is not surprising, therefore, that the 
response of disks is sensitive neither to the Reynolds number, nor to the amplitude of the 
initial disturbances. 

The mechanism by which shear flow nonlinearly breaks down into turbulence seems to 
be fairly-well understood, at least in its essentials: nonlinear perturbations are in essence 
linear disturbances that feed off of a marginally stable, finite amplitude two-dimensional per- 
turbation that is imposed on the shear flow (Zahn, Toomre, & Spiegel 1974; Bayly, Orszag, 
& Herbert 1988). If no such marginally stable, finite amplitude state exists — and for Keple- 
rian disks it does not — the mechanism breaks down. In fact, detailed numerical simulations 
show no tendency for Keplerian rotation to evolve toward self-sustained turbulence, even 
when driven vigorously well into the nonlinear regime by a sympathetic programmer. A 
revealing example is the case of thermally-driven convective turbulence. In the Boussinesq 
limit, this results in a very small but measurable angular momentum flux, with the wrong 
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sign to sustain turbulence (Stone & Balbus 1996). 1 Moreover, even codes with very different 
small scale dissipation properties converge when presented with the same initial set of finite 
amplitude disturbances (Hawley et al. 1999). Such convergence would hardly be expected 
if each code was anomalously dissipative in its own way. 

The current work is supportive of numerical simulation and the recent Ji laboratory 
experiment. The analytic perturbation solutions are in fact exact solutions to the dynamical 
equations. The nonlinear terms that are supposed to trigger turbulence instead vanish. This 
is problematic for the idea that these disturbances seed turbulence. More positively, it is 
an opportunity to test the accuracy of current codes. The solutions may be determined 
to any desired order of accuracy either by well-established WKB techniques, or by direct 
numerical integration of a simple ordinary differential equation. The wave forms involve time- 
dependent wavenumbers of arbitrarily high magnitudes, both leading and trailing. Despite 
the fact that the amplitude of the departure from Keplerian flow can be made arbitrarily high, 
the solution remains formally valid. Secondary instabilities feeding off this finite amplitude 
wave cannot be excluded a priori, but in fact the recent findings of Shen, Stone, & Gardiner 
(2006) suggest that such secondary instabilities make matters even more difficult. Very 
tightly-wound leading wave solutions (which would evolve toward the largest amplitudes) 
break down in the early stages of their evolution by what seems to be a Kelvin-Helmholtz 
instability. Far from becoming nonlinear, the wave dissipates before there is any appreciable 
growth. 



3. Preliminaries 

3.1. The Shearing Box 

The shearing box is a mathematical limit of the equations of motion in cylindrical coor- 
dinates that corresponds to retaining the local behavior of differentially rotating flow. The 
approach has a long history, dating at least as far back as the classical two-dimensional 
shearing sheet calculation of Goldreich & Lynden-Bell (1965). A more recent MHD descrip- 
tion can be find in the review article of Balbus & Hawley (1998). The idea is that one looks 
at the flow very near a fixed cylindrical radius R, and takes the limit that R — > oo. The 
angular velocity Q remains finite, and the azimuthal velocity formally increases without 



1 A negative Reynolds stress, which cannot be due to numerical dissipation, poses difficulties for the notion 
that Keplerian disks are formally turbulent, but at a level too small to be detected (Lesur & Longaretti 2005). 
Shear turbulence, whatever its magnitude, can be sustained only by a positive Reynolds stress. 
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bound. We are interested in departures from the Keplerian profile v^ — vk-, defined by 

u = v - v K e<p, (2) 

where v is the exact velocity, and is a unit vector in the azimuthal direction. Note that 
the u velocities are taken relative to the radially changing Keplerian profile, not relative to 
a frame that is rotating at a constant angular velocity. The unperturbed flow is thus u = 0, 
and | it | is always small compared with the large quantity v$. 

In a standard R, 0, z cylindrical coordinate system, the inviscid equations of motion in 
the local shearing box limit are (Balbus & Hawley 1998): 

(3) 

(4) 
(5) 
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(6) 

Our notation is as follows: p is the mass density of the disk, P is the gas pressure, and k 2 
is the epicyclic frequency associated with departures from circular orbits (eqn. 1). Formally, 
the local vertical gravitational force is —zQ 2 for a Keplerian potential. However it is often 
ignored in shearing box calculations (the "cylindrical approximation" ) and we shall generally 
do so in this paper. 

As noted, the equations apply to a local radial neighborhood around a large radius, R. 
The coefficients 2Q, k, 2 /2Q, and Q 2 are regarded as local constants. Finally, the equation for 
mass conservation in our constant density system is simply 

V-u = 0. (7) 

It is to be emphasized that these equations, while local, are fully nonlinear. 

It is also informative to write the shearing box equations in a frame that is rotating at 
a fixed angular velocity. Fix R to be the radius at the center of the box, and fix Q = Q(R). 
Let x be the radial distance from R, y = R(<f) — Qt), and z the vertical distance from the 
midplane. With 

w = v — Rile^, (8) 
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in their most general form the shearing box dynamical equations are 

- + w-Vj w R - 2Slw+ = - x — - --, (9) 
d \ ldP 

ai + w - v )*" + 2a »* = - f e? (10) 

l + w .v)w, = -zn>-~. (ii) 

at J p oz 

(We have retained the vertical force in equation [11].) Mass conservation is 

^ + V.( P w) = 0. (12) 

In this form, equations (9-12) are sometimes referred to as the Hill equations. 
The Hill equations have a simple scaling symmetry: 

w(r,t) <-> ew(r/e,t), (13) 

p(r,t) ~p(r/e,t), (14) 

P(r,t) <-> e 2 P(r/e,t). (15) 

In this form, with e < 1, we see that any solution of the Hill equations that involves 
very small length scales has a rescaled magnified counterpart with exactly the same time 
dependence. Any instability at small scales would also have to be unstable at larger scales. 



3.2. Shearing Coordinates 

We next introduce the coordinate transformation (Goldreich & Lynden-Bell 1965): 

(j)' = <f)-nt = (f)-t[n(R) + x(m/dR], (16) 

where we have expanded Q in a Taylor series about radius R (x is a small excursion from 
a fiducial value of R). The other coordinate transformations are identities, t' — t, Rl — R, 
z' = z. The transformed partial derivatives are 

d d d 
d d dVL d 

dR ~ dR> ~ t dR^ n ^ 
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dz dz r d(p dcf) 1 ' 
The explicit form of the equations of motion is 

du R 



du R „ „ 1/9 ,dn d \ „ 



(19) 



(20) 



where 



_i + u . Vl „ +_„„ = --_, (21) 

lF + "- v "' = -pS?' (22) 

"■ v ="«b-W + i* + v (23) 



Mass conservation is 



(9 9 \ 1 dus du z , ^ 



pB! dRd(j)'J Rdcj)' dz 
These equations form the basis of the analysis of the next section. 

4. Analysis 

4.1. Reduction to a Single ODE 

The system (20)-(24) has no spatial dependence, but does depend explicitly on time t'. 
This suggests that we seek a solution of the form 

fit') exp [i (k' R R' + m(f)'/R' + k z z')} , (25) 

where the function / will of course differ for each of the velocity components and the pressure. 
Without loss of generality, we may take m > 0. Notice that because of vanishing velocity 
divergence, the nonlinear (u • V)it terms vanish explicitly. Indeed, the u • V operator acting 
upon any perturbed modal plane wave vanishes. One can immediately see why nonlinear 
rotational hydro-instabilities are problematic: linear theory and nonlinear theory are the 
same for a finite amplitude plane wave. 

Although a superposition of these shearing waves would not be a solution of the nonlin- 
ear equations, in order to sustain turbulence such a superposition would have to somehow 
access the free energy source of differential rotation. But in general, the neglected nonlin- 
ear (u • 'V)u wave term leads not to a viable free energy source, but to mutual wave- wave 
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interactions. Absent a steady source term directly linked to the differential rotation, it is 
perhaps not surprising that these wave-wave cascades do not seem to be self-sustaining. 

In what follows, the fluid velocity u and pressure P variables refer to the time-dependent 
fourier / amplitude of equation (25). We define the unprimed time- dependent wavenumber 

k R = k R {t') = k' R -t'm^ (26) 

and shall henceforth drop the prime ' superscript from t' and other coordinates which are 
identical in both frames, and write d/dt for the time derivative of the amplitudes. The 
equations of motion are 

^ _ 2 n U4> = -ik R -, (27) 

du<i> k 2 ,mP ,_ x 

~i + ™ u « = ~'rj- (28 > 

du z , P , , 

= -ik z -, (29) 
dt p 

and the equation of mass conservation is 

u R k R + mu < p/R + k z u z = 0. (30) 

An equation for u R can be extracted from equations (26)-(30) after some straightforward, 
but tedious, algebra: 



d 2 u R 4k R m dVt du R 
~dt 2 k 2 R d\nR dt + 



k 2 o 2m 2 ( dtt 

t4« + 



k 2 k 2 R 2 \d\nR 



u R = 0, (31) 



where 



k 2 = k R (t) + ^ + k 2 . (32) 

This equation is in fact mathematically equivalent to equation (56) of Johnson & Gammie 
(2005a). These authors have shown that it may be solved in terms of hypergeometric func- 
tions. Here, we prefer to take a somewhat different tack, which avoids the introduction of 
rather complicated special functions. A pleasant simplification ensues if we introduce the 
variable 

U = k 2 u R . (33) 
Then, our differential equation takes on the compact Schrodinger form: 

d 2 ll k 2 

J + ¥ u - < 34 > 

and is amenable to a classical WKB analysis. 
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4-1.1. Explicit Solutions 

Equation (34) can be written in a convenient dimensionless form. Let 

k\ = k 2 + (m/R) 2 , k ± > 0, (35) 
and introduce a new time variable r: 

_ k R (t) 



k± 

Since dfl/dR < for the problems of interest, 

dr m dVt 



(36) 



, o > °- (37) 
dt k ± dR v 1 



With this substitution, our differential equation becomes 



dr 2 1 + r 2 

with 

P 



d 2 U B 2 

" ^ + T^-.U = 0, (38) 



k z K,/m 



(39) 



dtt/dR 

The leading order WKB solution of this equation is 

U ~ (1 + t 2 ) 1/4 exp [±ifi sinh" 1 r] , (40) 



or 



UR = u R (0) (1 + r 2 ) 1 exp [± i(3 sinh" 1 r] , (41) 

where ur(0) is defined to be ur at t — 0. Formally, this breaks down in the limit r ^> j3, 
but in practice very little changes. The differential equation for U may be easily solved to 
leading order in the overlapping asymptotic domain r ^> 1. For /3 > 1/4: 

u R = constant x r~ 3/2 exp [± - 1/A(3 2 ) 1/2 In r] , (42) 

while for < 1 /4: 

Mi? = constant x r^v 7 ^ )/2_ ( 43 ) 

With the precise leading order large r solution in hand, we see that equation (41) is actually 
valid in the domain lnr <^ 8/3. For (3 — 2, for example, this corresponds to r < 10 7 , 
considerably larger than the naive domain of applicability, r ^ f3 = 2. Even beyond the 
extended domain, however, the WKB solution differs from the actual large r solution only 
by a slowly varying phase. We shall henceforth adopt equation (41) as our formal solution 
for Ur. 
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4-1.2. Solution for the Pressure and Nonradial Velocities 



To complete the problem, we require expressions for u^, u z , and P. Unlike ur, these flow 
variables need not have a simple WKB form. Nevertheless, they may be obtained explicitly 
in terms of ur(t) and du^/dr. 



The key is to note that equation (38) may be written 

1 d 2 
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(3 2 dr 2 



which implies 



/X d 
u r (t) dr = -jp-^. [(l + t2 ) u r] + constant. 

If we combine equations (28), (29), and (30), there results 
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Writing this in terms of r and simplifying, 
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Finally, we integrate over r and use equation (45). We find 
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where the integration constant on the right has been set to zero to be consistent with the 
equation of motion (27). This determines u<j, in terms of ur and its first derivative: 
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Next, the axial velocity u z is obtained from V-u = 0: 
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The enthalpy P/p then follows from (27): 
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where it is understood that the real part of the equation is to be taken. This is convenient if 
u R is given as a WKB wave of the form (41), which is precisely how we shall use it. Choosing 



ur = Re 
we define the argument 



wr(0) exp[i(3 sinh (r)] 



ur(0) cos[/3sinh (r)] 



(1 + r 2 ) 3 / 4 

X = (3 sinh" 1 (t), 



(1 + r 2 ) 3 / 4 



and for Keplerian flow one obtains 
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(1 + r 2 ) 7 /4 



The velocity profiles are shown in Figure (1). 



(52) 
(53) 

(54) 

(55) 
(56) 



4.2. The Small (3 Limit 

The presence of the f3 term is a strongly stabilizing influence (it embodies epicyclic 
motion), and one therefore should examine the conditions under which it vanishes. In a 
Keplerian disk k = Q, and thus (3 1 must imply k z R m, which can correspond to one 
of two different physical situations. The first is when u z is a constant and all other variables 
vanish. These strictly vertical displacements are of little interest. The second possibility 
corresponds to degenerate incompressible two-dimensional disturbances in the disk plane. 
This problem has been considered by other authors (Chagelishvili et al. 2003; Johnson & 
Gammie 2005a). When (3 — > 0, we drop the second term in equation (38), and find the 
formal result 

Ci + C 2 r 

UR ~ I + t-2 ' ( } 

where C\ and C2 are constants of integration. This is consistent, of course, with equation 
(43) for large r. But in this limit equations (27), (28), and (30) may be directly combined 
into the first order equation 

M = &*L = 0t (58) 

GST GST 

so that C 2 must be zero. (It represents a spurious solution.) Chagelishvili et al. suggest 
that the growing phase of these solutions (t < 0) may lead to turbulence, but we have seen 



5 10 15 20 

T 

Fig. 1. — Radial (ur), azimuthal (u^), and axial (u z ) velocity components for disturbances 
in a (3 = 2 Keplerian disk. For reference, the (dashed line) radial velocity WKB solution lies 
almost directly on top of the exact solution ur, appearing slightly above the true solution 
near the minimum point. 
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that this solution is in fact exact. It cannot, and in simulations does not, lead to sustainable 
turbulence. 

All of this stands in stark contrast to the case of a vanishing specific angular momentum 
gradient (k 2 = 0), which does in fact lead to turbulence. Despite the apparently similar 
mathematics, the physical behavior of this (3 = solution is very different. Why? 

Because in the case of k 2 = 0, finite k = k z e z axisymmetric modes now are algebraically 
unstable. Without the steady march toward arbitrarily large radial wavenumbers that cuts- 
off nonaxisymmetric perturbations, axisymmetric ur perturbations continue to grow linearly 
with time. This solution is no longer a spurious branch of equation (34), it is a valid solution 
of the m = k = limit of equations (27)-(30). (Note that (3 is actually undefined in this 
singular limit, and r becomes a constant; one must work directly with t.) Perturbations 
in 11^ remain constant if m = 0, and grow linearly with time if the disturbance is even the 
slightest bit nonaxisymmetric. For a single mode, these are still exact nonlinear solutions, 
but an ensemble of modes in this case does lead to a turbulent cascade. The difference here 
is that the cascade can be fed by the axisymmetric instability, which features a consistently 
positive value for —URU^dQ/dR, and thus a robust source of free energy. Indeed, in the 
axisymmetric limit, equations (27) and (29) combine to give 

1 k du d , ., 

= mv^r > (59 > 



z 



so that outward angular momentum transport accompanies the release of free energy, growing 
linearly with time in the early stages. By contrast, non- axisymmetric solutions can maintain 
a source of free energy only for leading waves; this source becomes a sink when the waves 
become trailing. Unlike the uniform angular momentum turbulence, which is self-sustaining, 
fresh leading waves are always needed in a Keplerian disks. Without them, perturbations 
die. 



5. Simulations 

Numerical simulations have sharpened our understanding of the stability of Keplerian 
disks. These calculations are not infallible, of course: they are subject to the limitations that 
accompany any solution by numerical approximation. But the uncertainties due to numerical 
errors can generally be minimized by resolution studies, a careful choice of problem, and a 
detailed comparison with analytic solutions. This is as true for the hydrodynamic shearing 
box as for any other problem, although there are not many known analtyic solutions that 
can serve as tests. The shearing wave solutions presented in this paper constitute one such 
test problem. 
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Recently, the two-dimensional form of the shearing wave has been simulated by Johnson 
& Gammie (2005b), and by Shen, Stone & Gardiner (2006). These two studies used two 
different simulation codes: Johnson & Gammie (2005b) used the operator-split, nonconser- 
vative finite difference ZEUS algorithm (Stone & Norman 1992), and Shen et al. (2006) used 
the conservative Godunov Athena algorithm (Gardiner & Stone 2006). Here, we have carried 
out a few simple simulations using the Athena code to test the three-dimensional solutions 
developed in this paper. We have examined the evolution of a single linear shearing wave 
using the piecewise parabolic (third-order) representation of the underlying variables, and 
the Harten, Lax, van Leer with Contact (HLLC) flux option of the Athena code (Gardiner 
& Stone 2006). 

In general, a wave propagated by a finite difference algorithm will be subject to numerical 
truncation error, which is typically proportional to some power of the ratio of the grid 
resolution divided by the wavelength. The Athena code has a formal truncation error that is 
second order, and tests carried out by Shen et al. (2006) with two-dimensional shearing waves 
showed convergence to the analytic solution at a rate of iV~ 2 5 over a range from iV = 64 to 
256 grid zones. An interesting aspect of this test problem is that the radial wavenumber of a 
shearing wave evolves with time as the wave becomes increasingly tightly- wrapped (eq. 26). 
Thus, the effective resolution of any localized disturbance declines with time, regardless of 
the number of grid zones used. Shen et al. (2006) found that once a wave passed to fewer 
than 16 grid zones per wavelength, its amplitude decayed noticeably and smoothly as the 
resolution further diminished. Since the dissipation in Athena is strongly nonlinear a wave 
will be damped as its wavelength approached 2Ax, where Ax is a fiducial grid interval; Shen 
et al. observed no aliasing of trailing into leading waves. 

The specific problem we have examined is a Keplerian shearing box with q = —3/2, 
Q = 0.001, constant initial density p — 1, using both an isothermal equation of state with 
sound speed c s = 4.08f2, and an adiabatic equation of state with 7 = 5/3 and P = 10~ 6 . The 
computational domain has unit lengths, L R — — L z — 1. In local Cartesian coordinates 
dx = dR, dy = Rdcfi, the initial condition for the shearing wave is 

u R = A R sin[27r(n jR a; + n^y + n z z)} (60) 

where n« is the number of wavelengths within the domain: at t — 0, n R — 0, = 1, and n z = 
4, corresponding to /3 — 4. This initial condition is particularly convenient: and u z are zero 
at t — 0. The initial wave amplitude A R is set to 10~ 5 for the adiabatic simulation. We run 
the problem using 64 3 , 128 3 , and 256 3 grid zones. Figure (2) shows the evolution of the three 
velocity components for this mode in the 256 3 simulation overlaid on the analytic solution 
(dashed line). Some systematic deviation from the analytic solution becomes noticeable 
beyond n R = 16, which corresponds to 16 grid zones per radial wavelength for the 256 3 grid 
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zone resolution. The presence of small amplitude, high frequency sound waves can also be 
seen. 

The effect of resolution is illustrated in Figure (3) which plots the wave kinetic energy as 
a function of time for each of the three grid resolutions along with the analytic value (dashed 
line). As one would expect, as resolution increases the numerical solution maintains close 
adherence to the analytic value at higher n R values. Comparison of the fractional error for 
hr < 10 shows that between 128 3 and 256 3 zones the results are converging at the expected 
second order rate. The isothermal equation of state simulations are similar, although they 
deviate from the analytic value at an earlier time, suggesting that some of the difference is 
due to compressibility effects in addition to numerical sources of error. 

One reason that these shearing waves were of interest in previous two-dimensional stud- 
ies was that they are amplified during their leading phase. They reach a peak amplitude 
when r = 0; beyond that point they are trailing, and their amplitude declines. Numerical 
simulations offer a way to investigate the effects of compressibility and the stability of the 
waves as their peak velocity amplitudes approach the level of the background shear or the 
sound speed. Using 1024 2 grid zone two-dimensional simulations, Shen et al. (2006) exam- 
ined large amplitude waves and found that, for either leading or trailing waves, the waves 
become unstable when 

\k RU(t> \ > \dn/d\nR\. (61) 
When this limit is exceeded the waves rapidly break down. 

As a simple study of large amplitudes in three dimensions we start the (3 — 4 shearing 
wave mode at r = as above, but with a large amplitude, namely A R = 0.002. At this 
amplitude the stability criterion (61) will be violated when Su^ reaches its peak at n R ~ 2. 
This wave was evolved using 64 3 , 128 3 and 256 3 grid zones. Figure (4) shows the evolution of 
the kinetic energy in the velocity perturbations as a function of time. All three resolutions 
show a breakdown of the wave solution, although with increasing resolution the energy 
tracks the analytic solution for a slightly longer period of time before breaking down. The 
wave undergoes what appears to be a Kelvin- Helmholtz instability (Shen et al. 2006). The 
radial and azimuthal velocity amplitudes drop rapidly, leaving mainly vertical motions that 
have a finite ur but are nearly axisymmetric and ^-independent. This flow subsequently 
breaks down, again in a manner resembling a Kelvin-Helmholtz instability, and the wave 
kinetic energy rapid decays. Interestingly, the higher resolution simulation shows the most 
rapid decline in perturbed kinetic energy. Thus, while the high resolution simulations track 
the analytic solutions longer, they breakdown faster. This is perhaps a reflection of the 
fact that the highest resolved wavenumbers have the most rapid Kelvin-Helmholtz growth 
rates. The lower resolution simulations take longer for the secondary instability to affect 
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Fig. 2. — Amplitudes for radial (ur), azimuthal («</,), and axial (u z ) velocity components for 
a (3 = 4 wave, computed in a Keplerian shearing sheet using 256 3 grid zones. The dashed 
lines are the exact solutions derived from direct integration of the ODE. 
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Fig. 3. — Kinetic energy (as a fraction of the input value) in the velocity perturbations 
for a (3 = 4 Keplerian shearing sheet simulation. The labeled solid curves correspond to 
simulations using 64 3 , 128 3 and 256 3 grid zones. The dashed line is the analytic solution for 
the kinetic energy. 
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the vertical columns. As with the two-dimensional simulations, at late times the perturbed 
kinetic energy drops toward zero as does the (positive) Reynolds stress. There is no hint of 
sustained turbulence. 



6. Discussion and Conclusion 

We have argued that the growing behavior of the r < phase of local disk perturbations 
that has been put forth as a possible origin of sustained Keplerian turbulence will not, in 
fact, lead to such behavior. It is true that an arbitrarily large leading wavenumber will 
formally produce an arbitrarily large amplification factor during this phase of the evolution. 
The difficulty is that the linear and nonlinear behavior of the equations is the same for these 
incompressible plane wave solutions. The fact that the amplitudes become temporarily large 
does not change their wave character. Moreover, it appears that in the initial strongly leading 
configuration, the large shear in the wave is unstable, dissipating the disturbance at small 
amplitude before it has time to grow. 

True disk turbulence requires a tap into a permanent source of free energy, not large wave 
amplitudes. It is misleading to view the MRI and shear layer turbulence as a consequence 
of large amplitude perturbations. These flows are turbulent because there is a free energy 
path between the differential rotation and the perturbations, and it is already present at 
small amplitudes. The ensuing turbulent cascade is generally not be characterized by large 
amplitude disturbances. 

The outcome of numerical simulations form a coherent and interrelated pattern of re- 
sults. For example, Stone & Balbus (1996) studied the behavior of convective turbulence in 
Keplerian disks as a source of angular momentum transport. The correlation tensor (uru^), 
a measure of the angular momentum flux, was found to be tiny. Moreover, it had the wrong 
sign (radially inwards). Though not a direct investigation of Keplerian stability, these results 
make sense only if the underlying disk is stably stratified, and numerical diffusion is largely 
absent: reverse angular momentum transport is just what is expected in any Rayleigh-stable 
nonmagnetized disk that is driven externally (Balbus 2001). Detailed direct studies of Ke- 
plerian stability by Balbus et al. (1996), and Hawley et al. (1999) demonstrated nonlinear 
stability once more, even when the disks were strongly nonlinearly driven. Proponents of 
transient amplification as a route to turbulence must argue that the very large and tightly 
wound wavenumbers needed to enter into the nonlinear regime are inaccessible to current 
codes. But if this were the problem, one could simply start the calculation in the nonlinear 
regime. Flows that are nonlinearly unstable are in fact quite sensitive to starting ampli- 
tudes. When this is done for Keplerian disks, no transitions to sustained turbulence are 
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Fig. 4. — Kinetic energy (as a fraction of the input value) in the velocity perturbations for a 
j3 — 4 Keplerian shearing sheet simulation that began with a large initial amplitude. Shown 
are curves for 64 3 , 128 3 and 256 3 grid zones. The dashed line is the analytic solution for the 
kinetic energy. 
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found. Indeed, excellent convergence of all stable flows was found in the study by Hawley et 
al. (1999), even though the codes involved had completely different dissipative properties. 

The fact that the plane wave solution we have found is exact to all orders in its am- 
plitude suggests that coupling between the disturbance and the background is limited to 
simple wave action conservation. This is not a route to turbulence. A superposition of such 
plane wave solutions would exhibit group velocity dispersive spreading and nonlinear wave- 
wave interactions; the latter would impart power to smaller scales. Unless the fluid is itself 
externally driven, however, there is nothing to sustain the fluctuations against these losses. 

Because our solution is inviscid and three-dimensional, it provides an excellent diagnostic 
for probing the large dynamical range properties of numerical schemes. The solution is 
strictly exact only for constant density fluids, but the velocity components of an adiabatic 
fluid are accurately described by our solution. The large |r| behavior and peak amplitudes are 
stringent tests of a numerical scheme. The disagreement between the analytic and numerical 
large |r| behavior is particularly revealing, since it shows that the codes are able to uncover 
more complicated secondary features, even at very small scales. The most obvious effect 
in this work is compressibility, but the small scale Kelvin-Helmholtz instability reported by 
Shen et al. (2006) during the tightly-wrapped phase of the wave evolution (when \ur\ <C Iw^I) 
was unforeseen. A strongly leading shearing wave never even reaches the predicted analytic 
amplitude: its initial intrinsic shear is too unstable. 

Numerical, analytic, and laboratory studies appear to have reached a robust consensus: 
absent additional flow dynamics, Keplerian rotation is locally nonlinear stable. 
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